function [pval,W,df] = fsel_wald(X,Y,coef,r,set0,kmax)
% Forward selection for Wald test of H_0: beta(coef) = r(coef)
% Last edit: 3/28/17 CBH
% [p,W] = fsel_wald(X,Y,coef,r,set0,kmax)
% X - design matrix
% Y - outcome
% coef - variables whose coefficients are to be tested
% r - hypothesized value of coefficients
% set0 - set of initial regressors
% kmax - number of forward selection steps to take

p = size(X,2);
varind = (1:p)';

supp = set0;

W = zeros(kmax,1);
pval = zeros(kmax,1);
df = zeros(kmax,1);
for jj = 1:kmax
    varind_W = setdiff(varind,supp);
    
    Wjj = zeros(numel(varind_W),1);
    
    for kk = 1:numel(varind_W)
        if ~mod(kk,250)
            fprintf('Step: %d, Iteration: %d \n',jj,kk);
        end
        
        use = [supp;varind_W(kk)];
        z = X(:,use);
        b = z\Y;
        e = Y - z*b;
        [~,v] = hetero_se(z,e,inv(z'*z));
        B = zeros(p,1);
        V = zeros(p,p);
        B(use) = b;
        V(use,use) = v;
        
        testset = intersect(coef,use);
        Wjj(kk,1) = (B(testset)-r(testset))'*(V(testset,testset)\(B(testset)-r(testset)));
    end
    imin = find(Wjj == min(Wjj));
    imin = imin(1);
    
    supp = [supp ; varind_W(imin)]; %#ok<*AGROW>
        
    W(jj,1) = min(Wjj);
    pval(jj,1) = 1-chi2cdf(W(jj,1),numel(testset));
    df(jj,1) = numel(testset);
end
